clear all;
clc;
%%%%%%%%% Programming Language: Matlab 2019b       %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% (C) Andrew B. Abel and Stavros Panageas 2021 %%%%%%%%%%%%%%%%%%
%%%%%%%%% This code produces Figures 2,3 for the paper %%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% "An Analytic Framework for Interpreting %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% Investment Regressions in the Presence of Financial Constraints
%%%%%%%%% NOTE: The positioning of text labels in the figures may differ
%%%%%%%%% from the figures in the paper.
%%%%%%%%% This code also produces the entries for table 1.
%%%%%%%%% By setting report_only_table_3=1 and adjusting \alpha and \gamma
%%%%%%%%% the code 
%%%%%%%%% produces the entries for 
%%%%%%%%% table 3.
%%%%%%%%% By setting report_only_table_4=1 and adjusting the baseline parameters the code
%%%%%%%%% produces the entries for 
%%%%%%%%% table 3.


%%%%%%%%% Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
clc;

dt=0.01;

rho=0.07;%0.07
mu_H=0.3;%0.3
mu_L=0.1;%0.1
phi_L=-0.1;%-0.1
phi_H=0.2;%0.2
alpha=0.0; % Use alpha=0.2 for table 4. Adjust alpha to produce the entries of Table 3
gamma=0.0; % Adjust these values to produce Table 3 

report_only_table_3=0; % Set this to 1 to produce the entries of table 3.
report_only_table_4=0; % Set this to 1 to produce the entries of table 4.
rng(2525); % Fixes the random seed to allow reproducability of the results

theta=20; % This parameter is only used for the entries of table 4.
%%%% Derived parameter definitions and computation of x^ast %%%%%%%%%
            
eta_1=(rho+mu_H)/phi_L;
eta_2=(rho+mu_L)/phi_H;
G=mu_L/(rho+mu_L)*mu_H/(rho+mu_H);
Phi_L=1/eta_1;
Phi_H=1/eta_2;
% v_H_u=1/(1-G)*(Phi_H+mu_L/(rho+mu_L)*Phi_L);
% v_L_u=1/(1-G)*(Phi_L+mu_H/(rho+mu_H)*Phi_H);
% disp('Unconstrained cash flow coefficient')
% (v_H_u-v_L_u)/(phi_H-phi_L)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This section computes om_1 and om_2 and solves for x^ast
A=[eta_2 -mu_L/phi_H;-mu_H/phi_L eta_1];
om=eig(A);
om_1=om(2);
om_2=om(1);
x_m=max(1/(om_2-om_1)*log(om_1^2/om_2^2*(eta_2-om_2)/(eta_2-om_1)),0); %Compute the analytical formula of equation B.7 and use it as an initial guess
c_1=1/(om_2-om_1)*(om_2/om_1);
c_2=1/(om_1-om_2)*(om_1/om_2);
% The next equation corresponds to equation B.6 in the appendix:
mmf=@(y) c_1*(rho+mu_L)/mu_L*(1-om_1/eta_2)*exp(om_1*(-alpha*gamma-y))+c_2*(rho+mu_L)/mu_L*(1-om_2/eta_2)*exp(om_2*(-alpha*gamma-y)) - (1-gamma)*alpha;
disp('Critical Threshold: x^{\ast}')
x_m=fzero(mmf,x_m)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% Simulation start here %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
regime(1)=1;
nr_sim=1;
x(1)=0;
i=1;
cf(1)=phi_H;
while nr_sim<100 % This number controls how many repetitions to run. Increase this number for added precision
  flag=0; 
  
  if regime(i)==1 % Regime 1 is the "H" regime
        value(i)=c_1*exp(om_1*(x(i)-x_m))+c_2*exp(om_2*(x(i)-x_m));
        derivative(i)=om_1*c_1*exp(om_1*(x(i)-x_m))+c_2*om_2*exp(om_2*(x(i)-x_m));
        investment(i)=value(i)/derivative(i)-x(i); % As explained in the text, "investment" refers to marginal q.
        
    else
        value(i)=c_1*(rho+mu_L)/mu_L*(1-om_1/eta_2)*exp(om_1*(x(i)-x_m))+c_2*(rho+mu_L)/mu_L*(1-om_2/eta_2)*exp(om_2*(x(i)-x_m));
        derivative(i)=c_1*om_1*(rho+mu_L)/mu_L*(1-om_1/eta_2)*exp(om_1*(x(i)-x_m))+c_2*om_2*(rho+mu_L)/mu_L*(1-om_2/eta_2)*exp(om_2*(x(i)-x_m));
        investment(i)=value(i)/derivative(i)-x(i);
        
    end
    avq(i)=value(i)-x(i);
    
    
    
    if regime(i)==1 
        cf(i)=phi_H;
        x(i+1)=min(x_m, x(i)+phi_H*dt);
    else
        cf(i)=phi_L;
        x(i+1)=max(-gamma*alpha, x(i)+phi_L*dt);
    end
    
    if regime(i)==0 % If the firm goes bankrupt, restart the simulation with a new firm.
        if x(i+1)<=-gamma*alpha 
            nr_sim=nr_sim+1;
            regime(i+1)=1;
            flag=1;
            x(i+1)=0;
        end
    end
    
    
    
    
    if regime(i)==1 % Simulation of regime switching starts here
        if rand<1-exp(-mu_L*dt)
            regime(i+1)=0;
        else
            regime(i+1)=1;
        end
    else
        if flag==0
        if rand<1-exp(-mu_H*dt)
            regime(i+1)=1;
        else
            regime(i+1)=0;
        end
        end
    end
    i=i+1;
end
    
%%%% After simulation, run regressions on the simulated data.
% The output of these regressions are the numbers that show up in Table 3
% (after modifying alpha and gamma at the beginning of the file)
Y=investment';
X1=[avq' cf' ones(i-1,1)];
beta1=inv(X1'*X1)*X1'*Y;

if report_only_table_4==0
disp('This is the output for Table 3 and the first row of Table 1')
disp('Coefficient on Average q')
beta1(1)
disp('Coefficient on Cash Flow')
beta1(2)
else
disp('This is the output for Table 4 (Last two columns)')
disp('Coefficient on Average q')
beta1(1)/theta
disp('Coefficient on Cash Flow')
beta1(2)/theta
end





if report_only_table_3==0 && report_only_table_4==0
%%%%%%%%%%%% Plotting figures starts here %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ancillary computations to compute the weights of different observations %

m_x=length(x(x==x_m))/length(x);
[s0,xi0]=ksdensity(x(regime==0));
s0=s0/sum(s0);
[s1,xi1]=ksdensity(x(regime==1 & x<x_m));
s1=s1/sum(s1);
s=(regime.*interp1(xi1,s1,x)+(1-regime).*interp1(xi0,s0,x)).*(x<x_m)+m_x.*(x==x_m).*(regime==1);


%%%%%%%%%%%% UNCOMMENT THIS SECTION TO PRODUCE THE UNIVARIATE REGRESSION
%%%%%%%%%%%% RESULTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% X2=[avq' ones(i-1,1)];
% beta2=inv(X2'*X2)*X2'*Y;
% predicts2=max(X2*beta2,0);
% predicts2(find(predicts2==0))=NaN;
% figure
% title('Univariate regression - Figure 5')
% scatter(avq,investment,floor(s(1:end-1)*300)+1,[0 0 0]);
% hold on
% plot(avq,avq,'black')
% hold on
% plot(avq,predicts2,'red')


predicts3aa=max(X1*beta1,0).*(cf==phi_H)';
xvalsaa=avq.*(cf==phi_H);
predicts3aa(find(predicts3aa==0))=NaN;

predicts3bb=X1*beta1.*(cf==phi_L)';
xvalsbb=avq.*(cf==phi_L);

minq=min(avq);
maxq=max(avq);

xvalsa=0:0.001:maxq;
xvalsb=0:0.001:maxq;

predicts3a=beta1(3)+beta1(2)*phi_H+beta1(1)*xvalsa;
predicts3b=beta1(3)+beta1(2)*phi_L+beta1(1)*xvalsb;
predicts3a(find(predicts3a==0))=NaN;


%%%%%%%%%%%%%%%%%%%%%%%%%% Figure 2 (Note that text arrows might appear in the
%%%%%%%%%%%%%%%%%%%%%%%%%% wrong place)
%%%%%%%%%%%%%%%%%%%%%%%%%%
plot_additional=0;
figure
FS=14;
cics=0.03;
scatter(avq,investment,floor(s(1:end-1)*300)+1,[0 0 0])
hold on
text(min(avq),min(investment)+0.05,'O','FontSize',FS);
text(max(avq(regime==0)),max(investment(regime==0))+0.05,'B','FontSize',FS);
rectangle('Position',[min(avq)-cics/2,min(investment)-cics/2,cics,cics],'Curvature',[1 1],'FaceColor','red');
rectangle('Position',[max(avq(regime(1:end-1)==0))-cics/2,max(investment(regime(1:end-1)==0))-cics/2,cics,cics],'Curvature',[1 1],'FaceColor','red');
text(max(avq),max(investment)+0.05,'A','FontSize',FS);
text(min(avq(regime(1:end-1)==1))-0.05,min(investment(regime(1:end-1)==1))-0.05,'C','FontSize',FS);
rectangle('Position',[max(avq)-cics/2,max(investment)-cics/2,cics,cics],'Curvature',[1 1],'FaceColor','red');
rectangle('Position',[min(avq(regime(1:end-1)==1))-cics/2,min(investment(regime(1:end-1)==1))-cics/2,cics,cics],'Curvature',[1 1],'FaceColor','red');
plot(xvalsa,predicts3a,'--','Color','black')
hold on
plot(xvalsb,predicts3b,'--','Color','black')
hold on
plot(xvalsaa,predicts3aa,'.','Color','red')
hold on
plot(xvalsbb,predicts3bb,'.','Color','red')
xlabel('Average q','Interpreter','Latex','FontSize',FS)
ylabel('Marginal q','Interpreter','Latex','FontSize',FS)
box('off')
Xpos=0.5*(max(xvalsb)-min(xvalsa))*[1 1];
Ypos=0.55*[interp1(xvalsa,predicts3a,Xpos(1)),interp1(xvalsb,predicts3b,Xpos(1))];
annotation('doublearrow',[0.7 0.7],Ypos)
str={'Impact of cash flow', '$=(\phi^H-\phi^L)\beta_{\phi}$'};
textbo=annotation('textbox',[Xpos(1),Ypos(1),0.1,0.1],'String',str,'Interpreter','latex','EdgeColor','None');
textbo.Interpreter='latex';
textbo.FontSize=FS;
if plot_additional==1
    
    zz=get(gca,'XLim');
    hold on
    plot(linspace(zz(1),zz(2)*1.10,50),linspace(zz(1),zz(2)*1.10,50),'black')
    hold on
    X2z1=(linspace(zz(1),zz(2),50))';
    X2z2=[X2z1 ones(50,1)]*beta2;
    indsz=find(X2z2>zz(1));
    plot(X2z1(indsz),X2z2(indsz),'black')
    
    str={'45 degree line'};
    textbo=annotation('textarrow',[0.6,0.6],[0.6 0.5],'String',str,'Interpreter','latex');
    textbo.Interpreter='latex';
    textbo.FontSize=14;

    str={'Univariate regression on average $q$'};
    textbo=annotation('textarrow',[0.6,0.6],[0.6 0.5],'String',str,'Interpreter','latex');
    textbo.Interpreter='latex';
    textbo.FontSize=14;
    
    str={'Regression on average $q$ and cash flow'};
    textbo=annotation('textarrow',[0.6,0.6],[0.6 0.5],'String',str,'Interpreter','latex');
    textbo.Interpreter='latex';
    textbo.FontSize=14;
end


%%%%%%%%%%% TWO STAGE PROCEDURE STARTS HERE (Figure 3 and second row of Table 1) %%%%%%%%%%%% 


inds_regime_L=find(cf==phi_L);
Y=investment(inds_regime_L)';
X=[avq(inds_regime_L)' ones(length(inds_regime_L),1)];
beta4=inv(X'*X)*X'*Y;
predicts3c=beta4(2)+beta4(1)*xvalsb;
q_aa=max(avq(inds_regime_L));
q_mm=max(investment(inds_regime_L));
epsil=max(investment(inds_regime_L))-(beta4(2)+beta4(1)*max(avq(inds_regime_L)));

Th=(1-beta4(1))*maxq+epsil-q_mm+beta4(1)*q_aa; % Note that the max(q) and the max(marginal q) are equal in regime H
disp('Two-stage Procedure: Table 1 in the text')
[beta1(1) beta4(1); beta1(2) Th/(phi_H-phi_L)]'

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Figure 3 (Note that text arrows for the figure may appear at
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% the wrong places %%%%%%%%%%%%%%%%%%%%%%%%

figure
FS=14;
cics=0.03;
scatter(avq,investment,floor(s(1:end-1)*300)+1,[0 0 0])
hold on
text(min(avq),min(investment)+0.05,'O','FontSize',FS);
text(max(avq(regime==0)),max(investment(regime==0))+0.05,'B','FontSize',FS);
rectangle('Position',[min(avq)-cics/2,min(investment)-cics/2,cics,cics],'Curvature',[1 1],'FaceColor','red');
rectangle('Position',[max(avq(regime(1:end-1)==0))-cics/2,max(investment(regime(1:end-1)==0))-cics/2,cics,cics],'Curvature',[1 1],'FaceColor','red');
text(max(avq),max(investment)+0.05,'A','FontSize',FS);
rectangle('Position',[max(avq)-cics/2,max(investment)-cics/2,cics,cics],'Curvature',[1 1],'FaceColor','red');
plot(xvalsa,predicts3a,'--','Color','black')
hold on
plot(xvalsb,predicts3b,'--','Color','black')
hold on
plot(xvalsaa,predicts3aa,'.','Color','red')
hold on
plot(xvalsbb,predicts3bb,'.','Color','red')
xlabel('Average q','Interpreter','Latex','FontSize',FS)
ylabel('Marginal q','Interpreter','Latex','FontSize',FS)
box('off')
Xpos=0.5*(max(xvalsb)-min(xvalsa))*[1 1];
Ypos=0.55*[interp1(xvalsa,predicts3a,Xpos(1)),interp1(xvalsb,predicts3b,Xpos(1))];
annotation('doublearrow',[0.7 0.7],Ypos)
str={'$\Theta$'};
textbo=annotation('textbox',[Xpos(1),Ypos(1),0.1,0.1],'String',str,'Interpreter','latex','EdgeColor','None');
textbo.Interpreter='latex';
textbo.FontSize=FS;
hold on
plot(xvalsb,predicts3c,'Color','green');
text(max(xvalsb),max(predicts3c)+0.05,'D','FontSize',FS);
rectangle('Position',[max(xvalsb)-cics/2,max(predicts3c)-cics/2,cics,cics],'Curvature',[1 1],'FaceColor','red');
mr0=max(avq(regime==0));
line([mr0 mr0], [interp1(xvalsb,predicts3c,mr0) max(investment(regime==0))])
str={'$\zeta$'};
textbo=annotation('textarrow',[1.1*mr0 1.1*mr0],[0.3 0.4],'String',str,'Interpreter','latex');
textbo.Interpreter='latex';
textbo.FontSize=FS;
str={'Multiple regression (slope $\beta_q$)'};
textbo=annotation('textarrow',[0.6,0.6],[0.6 0.5],'String',str,'Interpreter','latex');
textbo.Interpreter='latex';
textbo.FontSize=12;
str={'Two-stage approximation (slope $\beta_q^{alt}$)'};
textbo=annotation('textarrow',[0.7,0.7],[0.2 0.3],'String',str,'Interpreter','latex');
textbo.Interpreter='latex';
textbo.FontSize=12;

end





